A comparison of axial and lateral PSF profiles of Field II against USTB's Fresnel simulator.

_by Alfonso Rodriguez-Molares alfonso.r.molares@ntnu.no, Ole Marius Hoel

Contents

Close old plots

close all;

Basic Constants

c0=1540;     % Speed of sound [m/s]
fs=100e6;    % Sampling frequency [Hz]
dt=1/fs;     % Sampling step [s]

field II initialisation

field_init(0);
set_field('c',c0);              % Speed of sound [m/s]
set_field('fs',fs);             % Sampling frequency [Hz]
set_field('use_rectangles',1);  % use rectangular elements
      *------------------------------------------------------------*
      *                                                            *
      *                      F I E L D   I I                       *
      *                                                            *
      *              Simulator for ultrasound systems              *
      *                                                            *
      *             Copyright by Joergen Arendt Jensen             *
      *    Version 3.30, April 5, 2021 (Matlab 2021a version)      *
      *                  Web-site: field-ii.dk                     *
      *                                                            *
      *     This is citationware. Note the terms and conditions    *
      *     for use on the web-site at:                            *
      *               field-ii.dk/?copyright.html                  *
      *  It is illegal to use this program, if the rules in the    *
      *  copyright statement is not followed.                      *
      *------------------------------------------------------------*
Warning:  Remember to set all pulses in apertures for the new sampling frequency

Transducer definition L9-4/38 Ultrasonix, 128-element linear array transducer

probe = uff.linear_array();
f0                      =5e6;             % Transducer center frequency [Hz]
lambda                  =c0/f0;           % Wavelength [m]
probe.element_height    =6e-3;            % Height of element [m]
probe.pitch             =0.3048e-3;       % probe.pitch [m]
kerf                    =0.035e-3;        % gap between elements [m]
probe.element_width     =probe.pitch-kerf;% Width of element [m]
lens_el                 =19e-3;           % position of the elevation focus
probe.N                 =128;             % Number of elements

Pulse definition

pulse = uff.pulse();
pulse.center_frequency = f0;
pulse.fractional_bandwidth = 0.6;             % probe bandwidth [1]
t0=(-1.0/pulse.fractional_bandwidth /f0): dt : (1.0/pulse.fractional_bandwidth /f0);
excitation=1;
impulse_response=gauspuls(t0, f0, pulse.fractional_bandwidth );
two_ways_ir= conv(conv(impulse_response,impulse_response),excitation)./norm(impulse_response).^2;
if mod(length(impulse_response),2)
    lag=(length(two_ways_ir)-1)/2;
else
    lag=(length(two_ways_ir))/2;
end

fig_handle=figure;
plot(((0:(length(two_ways_ir)-1))*dt -lag*dt)*1e6,two_ways_ir); hold on; grid on; axis tight
plot(((0:(length(two_ways_ir)-1))*dt -lag*dt)*1e6,abs(hilbert(two_ways_ir)),'r')
plot([0 0],[min(two_ways_ir) max(two_ways_ir)],'g');
legend('2-ways pulse','Envelope','Estimated lag');
title('2-ways impulse response Field II');
pulse.plot(fig_handle,'','--');

Aperture Objects

noSubAz=round(probe.element_width/(lambda/8));        % number of subelements in the azimuth direction
noSubEl=round(probe.element_height/(lambda/8));       % number of subelements in the elevation direction
Th = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]);
Rh = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]);

Excitation

xdc_excitation (Th, excitation);
xdc_impulse (Th, impulse_response);
xdc_baffle(Th, 0);
xdc_center_focus(Th,[0 0 0]);
xdc_impulse (Rh, impulse_response);
xdc_baffle(Rh, 0);
xdc_center_focus(Rh,[0 0 0]);

Phantom

pha=uff.phantom();
pha.sound_speed=1540;            % speed of sound [m/s]
pha.points=[0,  0, 20e-3, 1];    % point scatterer position [m]
fig_handle=pha.plot();
cropat=round(1.1*2*sqrt((max(pha.points(:,1))-min(probe.x))^2+max(pha.points(:,3))^2)/c0/dt);   % maximum time sample, samples after this will be dumped

Output data

t_out=0:dt:((cropat-1)*dt);                 % output time vector
STA=zeros(cropat,probe.N,probe.N);    % impulse response channel data

Compute STA signals using Field II

disp('Field II: Computing STA dataset');
wb = waitbar(0, 'Field II: Computing STA dataset');
for n=1:probe.N
    waitbar(n/probe.N, wb);

    % transmit aperture
    xdc_apodization(Th, 0, [zeros(1,n-1) 1 zeros(1,probe.N-n)]);
    xdc_focus_times(Th, 0, zeros(1,probe.N));

    % receive aperture
    xdc_apodization(Rh, 0, ones(1,probe.N));
    xdc_focus_times(Rh, 0, zeros(1,probe.N));

    % do calculation
    [v,t]=calc_scat_multi(Th, Rh, pha.points(1:3), pha.points(4));

    % build the dataset
    STA(1:size(v,1),:,n)=v;

    % Sequence generation
    seq(n)=uff.wave();
    seq(n).probe=probe;
    seq(n).source.xyz=[probe.x(n) probe.y(n) probe.z(n)];
    seq(n).sound_speed=c0;
    seq(n).delay = probe.r(n)/c0 - lag*dt + t; % t0 and center of pulse compensation
    seq(n).apodization = uff.apodization();
    seq(n).apodization.window=uff.window.sta;
end
close(wb);
Index exceeds the number of array elements (1).

Error in STAI_PSF (line 88)
disp('Field II: Computing STA dataset');

Channel Data

channel_data_field_ii = uff.channel_data();
channel_data_field_ii.sampling_frequency = fs;
channel_data_field_ii.sound_speed = c0;
channel_data_field_ii.initial_time = 0;
channel_data_field_ii.pulse = pulse;
channel_data_field_ii.probe = probe;
channel_data_field_ii.sequence = seq;
channel_data_field_ii.data = STA;

Scan

sca=uff.linear_scan('x_axis',linspace(-20e-3,20e-3,256).','z_axis', linspace(1e-3,40e-3,256).');

Pipeline

das = midprocess.das();
das.dimension = dimension.both;
das.channel_data=channel_data_field_ii;
das.scan=sca;

das.transmit_apodization.window=uff.window.hanning;
das.transmit_apodization.f_number=1.7;
%das.transmit_apodization.tilt=20*pi/180;
%das.transmit_apodization.plot()


das.receive_apodization.window=uff.window.hanning;
das.receive_apodization.f_number=1.7;
das.receive_apodization.tilt=20*pi/180;
%das.receive_apodization.probe = probe;
%das.receive_apodization.focus = sca;
%das.receive_apodization.plot()

b_data_field_ii =das.go();

b_data_field_ii.plot([],'Field II Simulation',70)

%das.receive_apodization.plot()